function V_hat = Vhat_NR(WTG, rho, CP, beta, lambda, LossData, RotSpd_act, Pow_act, PiPos2_act, V_hat_old)
% Define variables
R = WTG.R;
% Pnom = WTG.Pnom;
% Calculation of loss (if loss data used)
% [N1,M1]=size(LossData);
% Spd_ref=LossData(1,2:M1);   % [rpm], high speed side
% Load_ref=LossData(2:N1,1);  % [%], load/power relative to nominal
% Loss=LossData(2:N1,2:M1);   % [%], loss
% Loss1=interp2(Spd_ref,Load_ref,Loss,GenSpd_act,Pow_act/Pnom*100);
Loss1 = LossData; % to be outcommented if lossdata available

% check if pitch is in range of interpolation
if PiPos2_act > max(beta)
    PiPos2_act = 0.983 * max(beta);
    disp('Pitch out of range')
end
if PiPos2_act < min(beta)
    PiPos2_act = 0.9 * min(beta);
    disp('Pitch out of range')
end

% workaround to calculate the partial differential of Cp/lamda
dlambda = lambda(2) - lambda(1);
dCPdlambda = 1 / dlambda * diff(CP);
lambda_dCP = lambda(1:length(lambda) - 1) + 0.5 * dlambda; % this lambda vector must be used for interpolation in dCPdlambda

% Newton-Raphson iteration:
err = 1;      % error  
iteration = 0;
Vest = V_hat_old;    % estimated rotor wind initialisation (Vest=Vhat)

% Newton-Raphson root finding
while err > 0.01 && iteration < 100 % [m/s]
    lambda1 = RotSpd_act * R / Vest;
    % check if lambda is in range of interpolation
    if lambda1 > max(lambda_dCP)
        lambda1 = max(lambda_dCP);
        disp('Lambda out of range (max), V_hat')
    end;
    if lambda1 < min(lambda_dCP)
        lambda1 = min(lambda_dCP);
        disp('Lambda out of range (min), V_hat')
    end;
    CP1 = interp2(beta, lambda, CP, PiPos2_act, lambda1);
    dCPdlambda1 = interp2(beta, lambda_dCP, dCPdlambda, PiPos2_act, lambda1);
    % function f(v)
    % Gn(t)
    Fp = Pow_act / (Loss1 / 100) * 1e3 - 0.5 * rho * Vest^3 * pi * R^2 * CP1;
    % f'(v)
    % Gn'(t)
    dFp = 1/2 * rho * Vest^2 * pi * R^2 * (lambda1 * dCPdlambda1 - 3 * CP1);
    dFp = 1/2 * rho * pi * R^2 * Vest^2 * (lambda1 * dCPdlambda1 - 3 * CP1);
    % iteration form
    Vest_new = Vest - Fp / dFp;

    err = abs(Vest_new - Vest);
%     disp(err);
    Vest = Vest_new;

    iteration = iteration + 1;
end;
if iteration == 100
    disp(iteration);
end

% check gradient --> take out peaks (estimation errors)
if Vest > V_hat_old
    grad = Vest / V_hat_old;
else
    grad = V_hat_old / Vest;
end

if grad < 10
    V_hat = Vest;
else
    V_hat = V_hat_old;
    disp('Gradient exceeded')
end

if isnan(V_hat)
    V_hat = V_hat_old;   % if NaN for Vest set it to V_hat_old in order to run it
    disp('Interpolation of CP failed')
end;
